""" 
Author: H. Lin, Ph.D. 

First draft version: 2020;   Updated on 14th August 2025

R version 3.6.3 (2020-02-29) -- "Holding the Windsock"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
 
To cite R in publications use:   R Core Team (2020). R: A language and  environment for statistical computing. R   Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. 

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {R: A Language and Environment for Statistical Computing},
    author = {{R Core Team}},
    organization = {R Foundation for Statistical Computing},
    address = {Vienna, Austria},
    year = {2020},
    url = {https://www.R-project.org/},
  }


We have invested a lot of time and effort in creating R, please cite it when using it for data analysis. See also ‘citation("pkgname")’ for citing R packages.

"""







citation("e1071") # outputed below. Correct one, require "  " on "pack_name"

"""
> citation("e1071") 

  David Meyer, Evgenia Dimitriadou, Kurt Hornik,   Andreas Weingessel and Friedrich Leisch  (2019). e1071: Misc Functions of the
  Department of Statistics, Probability Theory   Group (Formerly: E1071), TU Wien. R package   version 1.7-3.   https://CRAN.R-project.org/package=e1071   

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {e1071: Misc Functions of the Department of Statistics, Probability
Theory Group (Formerly: E1071), TU Wien},
    author = {David Meyer and Evgenia Dimitriadou and Kurt Hornik and Andreas Weingessel and Friedrich Leisch},
    year = {2019},
    note = {R package version 1.7-3},
    url = {https://CRAN.R-project.org/package=e1071}, }              
"""





#  Data Preparation: Prepare and preprocess sequence data, vectors, machine learning labels, and splitting sets into the trainingset : cross-valiation set : testing set in the ratio of 6:2:2

#  Load the Human genome data into R environment
library(Biostrings) # Load the library. The function DNAString( ) of the library is required for format conversion purpose.

dna <- scan( # Load the sequence data into variable using scan( ) function
  
  file="sequence_file_path" , # Text file format for sequences 
  
  what = character(0) 
  
  ) 

class(dna) # Check the type of the variable. It was the character type

head(dna) # Display and check the head of the sequence data in the R environment

dna[1] #  Display and check the first line of the sequence data in the R environment

dna[2] #  Display and check the second line (data of the sequence in the R environment)
 
DNAString(dna[2]) # Convert the sequence text data to DNA string type data 

d3 <- DNAString(dna[2]); # Convert the sequence text data to DNA string type data and variable.  

class(d3) # Check the variable type. It was the "biostrings" type

d3[1];d3[2];d3[4]; substr(d3,1,1) # Display and check the bases in different numbered positions 







# Get the raw BRCA1 gene sequence from Human genome data
# Mark the beginning coordinate of BRCA1 gene in Gch38 version of human genome, which is 43044285
s1<-43044285; print(s1);  


df<-read.csv( 

  "Path_BRCA1_mutation_data_file.csv" ,  

  # From the disk file, Load BRCA1 mutation data into R environment, as a data.frame variable. The data include the mutatated bases' position number and the annotation of clinical signficances of those single point mutations. ( i.e., Whether a single base mutation is either benign or risky/pathogenic mutation. Such annotations are to be used as the label data in machine learning ) 

  header=T ) 
 

# Check the data structure of new variable  
str(df) # 1084 lines loaded into the data_frame type. Based on displayings, bases should be df[,4]; co-ordinates should be df[,3] 





mutants <- list() # Create an emtpy list for saving the mutated gene sequences 
 



for (i in 1:1084) {# A loop for creating gene BRCA1's mutant sequences. 1084 mutants in total.
  
  seq_op[df[i,3]-s1] <- as.character(df[i,4]); # replace one base each time. as.character( ) is required.

  mutants[i] <- seq_op; # save the mutant into a list var

  seq_op <- d3 # every time 's ending, recover it to wildtype/original standard BRCA1 seq is important!
};  



# To save/backup your data and variables in the R environment based on your needs. Could either do it via clicking graphical interface on Rstudio or use commands
# Command of Saving : save.image( )
save.image( # Saving data and variables in current R language computing environment/workspace
  
  file = "Path_your_data_file_name.RData", 
  
  version = 2.0,  safe = TRUE ) 



# To convert 1084 mutant sequences of BRCA1 gene into different vectors using different featurization methods (Such operation is also known as "embedding")

library(BioMedR) # Load the library required. Or use "require(BioMedR)" if data are reloaded from saved data file in the disk


v3mat <- matrix( # Build a matrix with specified and known dimension to store results.
  NA, 
  nrow = 1084, # Length/number of isntance
  ncol = 84 # Length/dimension of feature set
  )



for (i in 1:1084){ # a loop converting sequences data into k-mer vectors/values
  v3mat[i,] <- as.numeric(# Converting into numeric type variable
    extrDNAkmer(
      as.character(
        mutants[[i]]), # Specify the gene sequence ; Double "[[ ]]" must be used, for both variables
      k=3, # indicate 3-mer 
      upto=T # The argument "T" indicates both 1- , 2- and 3-mer vectors are generated.
      )
    ) # [[i]] must be used, for both variables
} ;  dim(v3mat) # cool! 1084 x 84 outputed


chk <- function(x){ # Build a function for checking the data properties of the vectors generated from above loop
  print(class(x)); # Display and check the class of variable
  print(head(x)) # Display and check the head content of the variable
  print(str(x)); # Display and check the data structure of variable
  print(summary(x)); # Display and check the data summary of variable 
  print(dim(x)); # Display and check the dimension of variable data
  print(length(x)) # Display and check the length of variable
  }; chk(v3mat); # Checking the data properties of the vectors generated from above loop




# To save/backup your data and variables in the R environment based on your needs. Could either do it via clicking graphical interface on Rstudio or use commands
# Command of Saving : save.image( )
save.image( # Saving data and variables in current R language computing environment/workspace  
  file = "Path_your_data_file_name.RData", 
  version = 2.0,  
  safe = TRUE 
  ) 


 




# Load the Support Vector Machine library "e1017" into R  environment. Or use "require(e1071)" to recover the library loading.
library(e1071) # require(e1071)

prime <- svm( # Specify inputs for the training of machine learning model (The Support Vector Machine here.)
  x=v3mat ,  # Specify the training set instance and feature set. Used full set of features here. 

  # It is possible to select / shorten the feature set to be used for machine learning. For example, "x = v3mat[ , 1:40]" tells the program to use only First to 40th columns as the machine learning training features. Similarly, "x = v3mat[ , X:Y]" tells the program to use only X-th to Y-th columns as the machine learning training features. 

  y = df[,2] , # Specify the label set for machine learning 
  
  cross = 10 , # Indicate the number of k for k-fold cross-valiation

  scale = T ,   # Indicate whether scaling is on or off.
  
  type = "nu-classification", #  Specify the purpose of machine learning, which type of classification or regression. Example options are  type="one-classification"), Or type="eps-regression" 
  
  kernel = "linear" # Specify the kernel of Support Vector Machine. Possible options are kernel = "polynomial",  kernel = "radial basis",   kernel = "radial",  kernel = "radial",    kernel = "sigmoid", etc. 
  
  ); # Note that, above arguments of cross, scale, type, kernerl can be ignored. When ignored, default arguments embedded in the function are to be used.

summary(prime) # Display the result summary of above 10-fold cross-validation of the primary machine learning, as shown below. But the performance was not good in the primary trial
"""
Call:
svm.default(x = v3mat, y = df[, 2], cross = 10)
Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 
Number of Support Vectors:  998
 ( 499 499 )
Number of Classes:  2 
Levels: 
 benign pathogenic

10-fold cross-validation on training data:
Total Accuracy: 53.96679 
Single Accuracies:
 52.77778 48.14815 54.12844 50.92593 54.12844 62.03704 63.88889 54.12844 50.92593 48.62385 
"""





    

